#'
#' calculate dispersion patterns of a population
#'
#'@param pop a population object
#'@param r vector of distance at which the metrics should be evaluated
#'@param metrics a vector of aggregation metrics should be calculated,see more information in details
#'@param correction a character to indicate which method is used to correct edge effect
#'@param approx whether to use samples of focus points to estimate the species dispersion
#'@param nample number of focus points used in approiximate estimation if approx is true.
#'
#'@details
#'This function is able to calculate six types of spatial characters of species distribution. specifically, 
#'the letters used in metrics mean:
#'
#'\tabular{ll}{
#' G(r) \tab the Nearest Neighbour Distance Function from an existed individual. 
#'      it is defined as a probability that the distance from a random chosen 
#'      event to its nearest neighbor is less than or equal to r. \cr
#' F(r) \tab the Nearest Neighbour Distance Function from an empty location. or empty space function,
#'      because it is a probability that the distance from a random chosen 
#'      empty location to its nearest neighbor is less than or equal to r. \cr
#' K(r) \tab the Ripley's K function. K-function is a second-moment measure as it is closely
#'     related to the second-order intensity of a stationary isotropic point process. 
#'     It captures the spatial (in)dependence between different regions of the point process.
#'      K-function is defined as: K(h) = l-1 E(# of other events within distance h of an arbitrary event). \cr
#' L(r) \tab the L function is just a simple sqrt(Kest/pi) transformation of the K function. \cr
#' g(r) \tab the pair correlation function. it defined as g(r) = K'(r)/ ( 2 * pi * r). \cr
#' O(r) \tab the O-ring function defined by lambda*g(r) \cr
#'}
#'
#'@note
#'the second and further high order nearest neighbor distances functions will be defined latter.
#'
#'@examples
#'set.seed(1)
#'
#'pop=rpopulation("sp1",n=5000,plotdim=c(100,100))
#'r=seq(0,5,length.out=20)
#'re=dispersion(pop,r,metrics=c("G","F"))
#'
#'data(BCI)
#'re=dispersion(BCI)
#'
#'bci=com_to_ppp(BCI)
#'Gest(bci,r=re$r)
#'
#'
#'#test the C++ version functions
#'
#'pop=rpopulation("sp1",n=50,plotdim=c(100,100))
#'re_r=dispersion(pop,r,metrics=c("G"))
#'

dispersion_C=function(data,r=NULL,metrics=c("G","F","K","L","g","O"),correction="border",approx=TRUE,nsample=1000){
  xy=as.matrix(data)
  plotdim=plotdim(data)
  xfocus=xy[,1]
  yfocus=xy[,2]
  #result is correct, but the speed is 105 times slower than the Gest(pop_to_ppp(pop),r,correction="none")
  #re=GFun(xy,plotdim[1],plotdim[2],xfocus,yfocus,r)
  #result is correct, but the speed is 135 times slower than the Kest(pop_to_ppp(pop),r=r,correction="none")
  #re=kFun(xy,plotdim[1],plotdim[2],xfocus,yfocus,r)
  #result is wrong, might caused by the normalization steps. Speed is 2 times faster than re<-pcf(pop_to_ppp(pop),r=r,correction="translate")
  re=gFun(xy,plotdim[1],plotdim[2],xfocus,yfocus,r)
  
  return(re)
}


dispersion=function(data,r=NULL,metrics=c("G","F","K","L","g","O"),correction="border",approx=TRUE,nsample=1000){
  plotdim=attr(data,"plotdim")
  lambda=total_abundance(data)/prod(plotdim)
  r_backup=r
  
  #define the length of r if it is not provided
  nr=50
  
  #if no approximate, just use all points
  if(!approx){
    nsample=total_abundance(data)
  }
  
  ##########################################
  #for G and F functions
  fx=numeric()
  fy=numeric()
  if("G" %in% metrics){
    #define the focus events
    fx=sample(data$x,nsample)
    fy=sample(data$y,nsample)
  }
  if("F" %in% metrics){
    #define the focus points
    fx=c(fx,runif(nsample,0,plotdim[1]))
    fy=c(fy,runif(nsample,0,plotdim[2]))
  }
  
  if(is.null(r)){
    rmax=rmax_default("G",plotdim,lambda)
    r=seq(0,rmax,length.out=nr)
  }else{
    rmax=max(r)
  }
  #get all neighborhood distances smaller than rmax
  alldist=frnn(data,fx,fy,rRange=c(0,rmax),type="circle",info="dist")
  #get the nearest neighborhood distances
  nndist=unlist(lapply(alldist,function(x) ifelse(length(x)!=0,min(x),rmax))) #get the nearest neighborhood distances
  
  re=data.frame(r=r)
  class(re)=c("dispersion",class(re))
  
  #calculate the empirical accumulate distribution
  if("G" %in% metrics){
    re$Gest=ecdf(nndist[1:nsample])(r)
    if(length(nndist)>nsample)
      nndist=nndist[(nsample+1):(2*nsample)]
  }
  if("F" %in% metrics){
    re$Fest=ecdf(nndist)(r)
  }
  
  if( !("K" %in% metrics) &  !("L" %in% metrics) & !("g" %in% metrics)){
   
    return(re)
  }
  
  #########################################
  #For K, L and g functions
  fx=numeric()
  fy=numeric()
  
  #define the focus events
  fx=sample(data$x,nsample)
  fy=sample(data$y,nsample)
  
  if(is.null(r_backup)){
    rmax=rmax_default("K",plotdim,lambda)
    r=seq(0,rmax,length.out=nr)
    re$r2=r
  }else{
    r=r_backup
    rmax=max(r)
  }
  
  #get all neighborhood distances smaller than rmax
  alldist=frnn(data,fx,fy,rRange=c(0,rmax),type="circle",info="dist")
  Kest=unlist(lapply(alldist,function(x) ecdf(x)(r)*length(x)))
  dim(Kest)=c(length(r),nsample)
  #the K function
  Kest=apply(Kest,1,mean)/lambda
    
  if("K" %in% metrics)
    re$Kest=Kest
  
  if("L" %in% metrics)
    re$Lest=sqrt(Kest/pi)

  
  if("g" %in% metrics | "O" %in% metrics){
    
    z <- Kest/(pi*r^2)
    z[!is.finite(z)] <- 1
    ss <- smooth.spline(x = r, y = z)
    dz <- predict(ss, r, deriv = 1)$y
    
    if("g" %in% metrics){
      re$gest <- (r/2) * dz + z
      if("O" %in% metrics){
        re$O=re$gest*lambda
      }
    }else{
      re$O=((r/2) * dz + z)*lambda
    }
    
  }
  
 

  return(re)
}


rmax_default=function (fun = "K", plotdim, lambda) 
{
  #verifyclass(W, "owin")
  switch(fun, K = {
    ripley <- min(plotdim)/4
    rlarge <- if (!missing(lambda)) sqrt(1000/(pi * lambda)) else Inf
    rmax <- min(rlarge, ripley)
  }, F = , G = , J = {
    rdiam <- sqrt(plotdim[1]^2+plotdim[2]^2)/2
    rlarge <- if (!missing(lambda)) sqrt(log(1e+05)/(pi * 
                                                       lambda)) else Inf
    rmax <- min(rlarge, rdiam)
  }, stop(paste("Unrecognised function type", sQuote(fun))))
  return(rmax)
}

